library("tidyverse")

# process plot data
d <- read.csv("output/summary.csv")
d$date <- as.Date(d$date)

### 1. Plot multianual timeseries
png(paste0("output/pics/chla-dynamics-summary.png"), res=180, width = 1200, height=850)
plot(d$date, d$df.chl_mean, type="l", main = "Чограйское водохранилище", ylim=c(0,25), lwd=2, xlab="Вектор времени", ylab="Хлорофилл-а, мкг/л")
lines(d$date, d$df.chl_min, lty=2, col="blue", lwd = 0.5)
lines(d$date, d$df.chl_max, lty=2, col="red", lwd = 0.5)
legend("topleft", legend = c("min", "average", "max"), fill=c("blue", "black", "red"))
dev.off()
#axis(1, d$date, format(d$date, "%Y-%m-%d"))

# calc grouped data by month by all years
md <- data.frame(m = d$df.month, v = d$df.chl_mean) %>%
  group_by(m) %>%
  summarise(mean = mean(v), n = n())

d21 <- d[d$df.year == 2021,]
md21 <- data.frame(m = d21$df.month, v = d21$df.chl_mean) %>%
  group_by(m) %>%
  summarise(mean = mean(v), n = n())

d22 <- d[d$df.year == 2022,]
md22 <- data.frame(m = d22$df.month, v = d22$df.chl_mean) %>%
  group_by(m) %>%
  summarise(mean = mean(v), n = n())

d23 <- d[d$df.year == 2023,]
md23 <- data.frame(m = d23$df.month, v = d23$df.chl_mean) %>%
  group_by(m) %>%
  summarise(mean = mean(v), n = n())

### 2. Plot seasonal density
png(paste0("output/pics/chla-seasonality.png"), res=180, width = 1200, height=850)
xlab = 1:12
plot(d$df.month, d$df.chl_mean, main="Сезонная динамика хлорофилла-а", xlab="№ месяца", ylab="Концентрация хлорофилла а, мкг/л", xlim=c(1,12), xaxt="n", ylim=c(0,12), cex=0.4)
axis(1, at=xlab)
lines(md$m, md$mean, lwd=3, col="black", lty=1)
lines(md21$m, md21$mean, col="darkgreen", lty=2)
lines(md22$m, md22$mean, col="darkblue", lty=2)
lines(md23$m, md23$mean, col="darkred", lty=2)
legend("topleft", legend = c("Среднее 2021-2023", "2021", "2022", "2023"), lty = c(1,2,2,2), col = c("black", "darkgreen", "darkblue", "darkred"))

dev.off()


### 3. Convert to pcb (phytoplankton biomass) by Минеева и др., 2014
# eq: pcb = 0,225 * chl(a)
d$pcb <- 0.225 * d$df.chl_mean
pcb <- data.frame(m = d$df.month, v = d$pcb) %>%
  group_by(m) %>%
  summarise(mean = mean(v), n = n())
write.csv(pcb, "output/pcb.csv", row.names = F)
